Edit on GitHub

McmcHermes.jl

A documentation for the McmcHermes package.

McmcHermes provides a simple but efficient way to generate Markov Chain Monte-Carlo algorithms in order to sample a probability density distribution.

Overview

The major functions in this module are:

runmcmc: run multiple chains with a specific number of walkers. get\flat_chain: get the stored chain of MCMC samples. get_gelman_rubin: get the Gelman Rubin convergence diagnostic of the chains.

Note

This guide assumes that you already have define your likelihood, prior and the logarithm of the posterior probability as in the example below.

Pkg Registry

using Pkg
Pkg.add("McmcHermes")

Example

First, let's generate some data:

using Distributions, Plots, LaTeXStrings, DataFrames, ProgressMeter

mu, sigma = 10, 2
l_b, u_b = 0, 20
d = Truncated(Normal(mu, sigma), l_b, u_b)
N = 1000
data = rand(d, N)

histogram(data, legend=false, size=(300,300), xlabel="data", show=true)

data

In order to sample the posterior probability distribution, it is necessary to define the likelihood, prior and logarithm of the posterior probability.

function log_likelihood(X::Vector, parameters::Vector)
    mu, sigma = parameters[1], parameters[2]
    y = 1 ./ (sqrt(2 * pi) .* sigma) .* exp.( -0.5 * ((X .- mu)./sigma).^2 )
    return sum(log.(y))
end

function log_prior(parameters::Vector)
    mu, sigma = parameters[1], parameters[2]
    if 5.0 < mu < 15.0 && 0.0 < sigma < 4.0
        return 0.0
    end
    return -Inf
end

function log_probability(X::Vector, parameters::Vector)
    lp = log_prior(parameters)
    if !isfinite(lp)
        return -Inf
    end
    return lp + log_likelihood(X, parameters)
end
log_probability (generic function with 1 method)

Call the McmcHermes package and define the number of walkers, iterations, dimension of the parameter space and the initial guess.

using McmcHermes

mu, sigma = 10, 2
initparams = Vector{Float64}([mu, sigma])

n_iter, n_walkers = 1000, 200
n_dim, a = 2, 0.01

chain_tests = run_mcmc(log_probability, data, initparams, n_iter, n_walkers, n_dim, a=a)
println(size(chain_tests))

Running chains...   1%|▍                                 |  ETA: 0:01:13
Running chains...   6%|██                                |  ETA: 0:00:43
Running chains...   6%|██▎                               |  ETA: 0:00:41
Running chains...   7%|██▍                               |  ETA: 0:00:40
Running chains...   8%|██▌                               |  ETA: 0:00:39
Running chains...   8%|██▊                               |  ETA: 0:00:38
Running chains...   8%|██▉                               |  ETA: 0:00:37
Running chains...   9%|███                               |  ETA: 0:00:36
Running chains...  10%|███▎                              |  ETA: 0:00:36
Running chains...  10%|███▍                              |  ETA: 0:00:35
Running chains...  10%|███▋                              |  ETA: 0:00:34
Running chains...  11%|███▊                              |  ETA: 0:00:34
Running chains...  12%|███▉                              |  ETA: 0:00:33
Running chains...  12%|████▏                             |  ETA: 0:00:33
Running chains...  12%|████▎                             |  ETA: 0:00:32
Running chains...  13%|████▍                             |  ETA: 0:00:32
Running chains...  14%|████▋                             |  ETA: 0:00:31
Running chains...  14%|████▊                             |  ETA: 0:00:31
Running chains...  14%|████▉                             |  ETA: 0:00:31
Running chains...  15%|█████▏                            |  ETA: 0:00:30
Running chains...  16%|█████▎                            |  ETA: 0:00:30
Running chains...  16%|█████▌                            |  ETA: 0:00:29
Running chains...  16%|█████▋                            |  ETA: 0:00:29
Running chains...  17%|█████▊                            |  ETA: 0:00:29
Running chains...  18%|██████                            |  ETA: 0:00:29
Running chains...  18%|██████▏                           |  ETA: 0:00:28
Running chains...  18%|██████▎                           |  ETA: 0:00:28
Running chains...  19%|██████▌                           |  ETA: 0:00:28
Running chains...  20%|██████▋                           |  ETA: 0:00:27
Running chains...  20%|██████▊                           |  ETA: 0:00:27
Running chains...  20%|███████                           |  ETA: 0:00:27
Running chains...  21%|███████▏                          |  ETA: 0:00:27
Running chains...  22%|███████▎                          |  ETA: 0:00:26
Running chains...  22%|███████▌                          |  ETA: 0:00:26
Running chains...  22%|███████▋                          |  ETA: 0:00:26
Running chains...  23%|███████▉                          |  ETA: 0:00:26
Running chains...  24%|████████                          |  ETA: 0:00:25
Running chains...  24%|████████▏                         |  ETA: 0:00:25
Running chains...  24%|████████▍                         |  ETA: 0:00:25
Running chains...  25%|████████▌                         |  ETA: 0:00:25
Running chains...  26%|████████▋                         |  ETA: 0:00:24
Running chains...  26%|████████▉                         |  ETA: 0:00:24
Running chains...  26%|█████████                         |  ETA: 0:00:24
Running chains...  27%|█████████▏                        |  ETA: 0:00:24
Running chains...  28%|█████████▍                        |  ETA: 0:00:24
Running chains...  28%|█████████▌                        |  ETA: 0:00:23
Running chains...  28%|█████████▊                        |  ETA: 0:00:23
Running chains...  29%|█████████▉                        |  ETA: 0:00:23
Running chains...  30%|██████████                        |  ETA: 0:00:23
Running chains...  30%|██████████▎                       |  ETA: 0:00:23
Running chains...  30%|██████████▍                       |  ETA: 0:00:22
Running chains...  31%|██████████▌                       |  ETA: 0:00:22
Running chains...  32%|██████████▊                       |  ETA: 0:00:22
Running chains...  32%|██████████▉                       |  ETA: 0:00:22
Running chains...  32%|███████████                       |  ETA: 0:00:22
Running chains...  33%|███████████▎                      |  ETA: 0:00:21
Running chains...  34%|███████████▍                      |  ETA: 0:00:21
Running chains...  34%|███████████▌                      |  ETA: 0:00:21
Running chains...  34%|███████████▊                      |  ETA: 0:00:21
Running chains...  35%|███████████▉                      |  ETA: 0:00:21
Running chains...  36%|████████████▏                     |  ETA: 0:00:20
Running chains...  36%|████████████▎                     |  ETA: 0:00:20
Running chains...  36%|████████████▍                     |  ETA: 0:00:20
Running chains...  37%|████████████▋                     |  ETA: 0:00:20
Running chains...  38%|████████████▊                     |  ETA: 0:00:20
Running chains...  38%|████████████▉                     |  ETA: 0:00:20
Running chains...  38%|█████████████▏                    |  ETA: 0:00:19
Running chains...  39%|█████████████▎                    |  ETA: 0:00:19
Running chains...  40%|█████████████▍                    |  ETA: 0:00:19
Running chains...  40%|█████████████▋                    |  ETA: 0:00:19
Running chains...  40%|█████████████▊                    |  ETA: 0:00:19
Running chains...  41%|██████████████                    |  ETA: 0:00:18
Running chains...  42%|██████████████▏                   |  ETA: 0:00:18
Running chains...  42%|██████████████▎                   |  ETA: 0:00:18
Running chains...  42%|██████████████▌                   |  ETA: 0:00:18
Running chains...  43%|██████████████▋                   |  ETA: 0:00:18
Running chains...  44%|██████████████▊                   |  ETA: 0:00:18
Running chains...  44%|███████████████                   |  ETA: 0:00:17
Running chains...  44%|███████████████▏                  |  ETA: 0:00:17
Running chains...  45%|███████████████▎                  |  ETA: 0:00:17
Running chains...  46%|███████████████▌                  |  ETA: 0:00:17
Running chains...  46%|███████████████▋                  |  ETA: 0:00:17
Running chains...  46%|███████████████▊                  |  ETA: 0:00:17
Running chains...  47%|████████████████                  |  ETA: 0:00:16
Running chains...  48%|████████████████▏                 |  ETA: 0:00:16
Running chains...  48%|████████████████▍                 |  ETA: 0:00:16
Running chains...  48%|████████████████▌                 |  ETA: 0:00:16
Running chains...  49%|████████████████▋                 |  ETA: 0:00:16
Running chains...  50%|████████████████▉                 |  ETA: 0:00:16
Running chains...  50%|█████████████████                 |  ETA: 0:00:15
Running chains...  50%|█████████████████▏                |  ETA: 0:00:15
Running chains...  51%|█████████████████▍                |  ETA: 0:00:15
Running chains...  52%|█████████████████▌                |  ETA: 0:00:15
Running chains...  52%|█████████████████▋                |  ETA: 0:00:15
Running chains...  52%|█████████████████▉                |  ETA: 0:00:15
Running chains...  53%|██████████████████                |  ETA: 0:00:14
Running chains...  54%|██████████████████▎               |  ETA: 0:00:14
Running chains...  54%|██████████████████▍               |  ETA: 0:00:14
Running chains...  54%|██████████████████▌               |  ETA: 0:00:14
Running chains...  55%|██████████████████▊               |  ETA: 0:00:14
Running chains...  56%|██████████████████▉               |  ETA: 0:00:14
Running chains...  56%|███████████████████               |  ETA: 0:00:13
Running chains...  56%|███████████████████▎              |  ETA: 0:00:13
Running chains...  57%|███████████████████▍              |  ETA: 0:00:13
Running chains...  58%|███████████████████▌              |  ETA: 0:00:13
Running chains...  58%|███████████████████▊              |  ETA: 0:00:13
Running chains...  58%|███████████████████▉              |  ETA: 0:00:13
Running chains...  59%|████████████████████              |  ETA: 0:00:12
Running chains...  60%|████████████████████▎             |  ETA: 0:00:12
Running chains...  60%|████████████████████▍             |  ETA: 0:00:12
Running chains...  60%|████████████████████▋             |  ETA: 0:00:12
Running chains...  61%|████████████████████▊             |  ETA: 0:00:12
Running chains...  62%|████████████████████▉             |  ETA: 0:00:12
Running chains...  62%|█████████████████████▏            |  ETA: 0:00:12
Running chains...  62%|█████████████████████▎            |  ETA: 0:00:11
Running chains...  63%|█████████████████████▍            |  ETA: 0:00:11
Running chains...  64%|█████████████████████▋            |  ETA: 0:00:11
Running chains...  64%|█████████████████████▊            |  ETA: 0:00:11
Running chains...  64%|█████████████████████▉            |  ETA: 0:00:11
Running chains...  65%|██████████████████████▏           |  ETA: 0:00:11
Running chains...  66%|██████████████████████▎           |  ETA: 0:00:10
Running chains...  66%|██████████████████████▌           |  ETA: 0:00:10
Running chains...  66%|██████████████████████▋           |  ETA: 0:00:10
Running chains...  67%|██████████████████████▊           |  ETA: 0:00:10
Running chains...  68%|███████████████████████           |  ETA: 0:00:10
Running chains...  68%|███████████████████████▏          |  ETA: 0:00:10
Running chains...  68%|███████████████████████▎          |  ETA: 0:00:09
Running chains...  69%|███████████████████████▌          |  ETA: 0:00:09
Running chains...  70%|███████████████████████▋          |  ETA: 0:00:09
Running chains...  70%|███████████████████████▊          |  ETA: 0:00:09
Running chains...  70%|████████████████████████          |  ETA: 0:00:09
Running chains...  71%|████████████████████████▏         |  ETA: 0:00:09
Running chains...  72%|████████████████████████▎         |  ETA: 0:00:09
Running chains...  72%|████████████████████████▌         |  ETA: 0:00:08
Running chains...  72%|████████████████████████▋         |  ETA: 0:00:08
Running chains...  73%|████████████████████████▉         |  ETA: 0:00:08
Running chains...  74%|█████████████████████████         |  ETA: 0:00:08
Running chains...  74%|█████████████████████████▏        |  ETA: 0:00:08
Running chains...  74%|█████████████████████████▍        |  ETA: 0:00:08
Running chains...  75%|█████████████████████████▌        |  ETA: 0:00:08
Running chains...  76%|█████████████████████████▋        |  ETA: 0:00:07
Running chains...  76%|█████████████████████████▉        |  ETA: 0:00:07
Running chains...  76%|██████████████████████████        |  ETA: 0:00:07
Running chains...  77%|██████████████████████████▏       |  ETA: 0:00:07
Running chains...  78%|██████████████████████████▍       |  ETA: 0:00:07
Running chains...  78%|██████████████████████████▌       |  ETA: 0:00:07
Running chains...  78%|██████████████████████████▊       |  ETA: 0:00:06
Running chains...  79%|██████████████████████████▉       |  ETA: 0:00:06
Running chains...  80%|███████████████████████████       |  ETA: 0:00:06
Running chains...  80%|███████████████████████████▎      |  ETA: 0:00:06
Running chains...  80%|███████████████████████████▍      |  ETA: 0:00:06
Running chains...  81%|███████████████████████████▌      |  ETA: 0:00:06
Running chains...  82%|███████████████████████████▊      |  ETA: 0:00:06
Running chains...  82%|███████████████████████████▉      |  ETA: 0:00:05
Running chains...  82%|████████████████████████████      |  ETA: 0:00:05
Running chains...  83%|████████████████████████████▎     |  ETA: 0:00:05
Running chains...  84%|████████████████████████████▍     |  ETA: 0:00:05
Running chains...  84%|████████████████████████████▌     |  ETA: 0:00:05
Running chains...  84%|████████████████████████████▊     |  ETA: 0:00:05
Running chains...  85%|████████████████████████████▉     |  ETA: 0:00:04
Running chains...  86%|█████████████████████████████▏    |  ETA: 0:00:04
Running chains...  86%|█████████████████████████████▎    |  ETA: 0:00:04
Running chains...  86%|█████████████████████████████▍    |  ETA: 0:00:04
Running chains...  87%|█████████████████████████████▋    |  ETA: 0:00:04
Running chains...  88%|█████████████████████████████▊    |  ETA: 0:00:04
Running chains...  88%|█████████████████████████████▉    |  ETA: 0:00:04
Running chains...  88%|██████████████████████████████▏   |  ETA: 0:00:03
Running chains...  89%|██████████████████████████████▎   |  ETA: 0:00:03
Running chains...  90%|██████████████████████████████▍   |  ETA: 0:00:03
Running chains...  90%|██████████████████████████████▋   |  ETA: 0:00:03
Running chains...  90%|██████████████████████████████▊   |  ETA: 0:00:03
Running chains...  91%|███████████████████████████████   |  ETA: 0:00:03
Running chains...  92%|███████████████████████████████▏  |  ETA: 0:00:03
Running chains...  92%|███████████████████████████████▎  |  ETA: 0:00:02
Running chains...  92%|███████████████████████████████▌  |  ETA: 0:00:02
Running chains...  93%|███████████████████████████████▋  |  ETA: 0:00:02
Running chains...  94%|███████████████████████████████▊  |  ETA: 0:00:02
Running chains...  94%|████████████████████████████████  |  ETA: 0:00:02
Running chains...  94%|████████████████████████████████▏ |  ETA: 0:00:02
Running chains...  95%|████████████████████████████████▎ |  ETA: 0:00:01
Running chains...  96%|████████████████████████████████▌ |  ETA: 0:00:01
Running chains...  96%|████████████████████████████████▋ |  ETA: 0:00:01
Running chains...  96%|████████████████████████████████▊ |  ETA: 0:00:01
Running chains...  97%|█████████████████████████████████ |  ETA: 0:00:01
Running chains...  98%|█████████████████████████████████▏|  ETA: 0:00:01
Running chains...  98%|█████████████████████████████████▍|  ETA: 0:00:01
Running chains...  98%|█████████████████████████████████▌|  ETA: 0:00:00
Running chains...  99%|█████████████████████████████████▋|  ETA: 0:00:00
Running chains... 100%|█████████████████████████████████▉|  ETA: 0:00:00
Running chains... 100%|██████████████████████████████████| Time: 0:00:29
(1000, 200, 2)

Gelman-Rubin's diagnostic can be obtained from the chains calling the get_gelman_rubin method.

println("Gelman Rubin Diagnostic: ", get_gelman_rubin(chain_tests))
Gelman Rubin Diagnostic: 1.057201558944357

Finally, plot the chains.

labels = Tuple([L"\mu", L"\sigma"])
x = 1:size(chain_tests)[1]
p = []
for ind in 1:n_dim
    push!(p, plot(x, [chain_tests[:,i,ind] for i in 1:size(chain_tests)[2]], legend=false,
    lc=:black, lw=1, ylabel=labels[ind], alpha = 0.1, xticks=false))
end

plot(p[1], p[2], layout = (2,1))
plot!(size=(600,200), xlims = (0, size(chain_tests)[1]), show=true)

chains

flat_chains = get_flat_chain(chain_tests, burn_in=100, thin=10)

flat = DataFrame(flat_chains, :auto)
colnames = ["mu", "sigma"]
flat = rename!(flat, Symbol.(colnames))

using PairPlots, CairoMakie
pairplot(flat)

corner

Develop by Steven Alfonso.